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Formation of young boxy/peanut bulges in ringed barred galaxies 
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ABSTRACT 



Aims. We investigate whether the formation mechanism of boxy and peanut-shaped (B/PS) bulges could depend on the gas content of 
the galaxy. 

Methods. We have performed N-body simulations with and without a gaseous component. In the second case star formation/feedback 
recipes have also been implemented to create new stellar populations. 

Results. As in many previous studies, in our N-body collisionless simulation, the B/PS is due to the classical break in the z mirror 
symmetry lasting roughly 200 Myr. When a gaseous component and star formation recipes are added to the simulation, the bulge- 
growing mechanism is quite different. The young stellar population that is born in the thin gaseous disc rapidly populates vertical 
resonant orbits triggered by the combined effects of the linear horizontal and vertical ILRs. This leads to a B/PS bulge mainly made 
of stellar material younger than the surrounding population. 

The non-linear analysis of the orbital structure shows that the main orbit family responsible for the B/PS is not the same in the two 
cases. The 2:2:1 orbits prevail in the collisionless simulation whereas additional asymmetrical families contribute to the B/PS if a 
dissipative component is present and can form new stars. We found that 2:3:1 and 2:5:1 orbits trap a significant fraction of the mass. 
A flat ringed discy stellar component also appears simultaneously with the thickening of the young population. It is due to the star 
formation in a nuclear gaseous disc located in the central kpc, inside the ILR, and accumulated there by the torques exerted by the 
large-scale bar. Remarkably, it remains flat throughout the simulation although it develops a nuclear bar, leading to a double-barred 
galaxy. 

Conclusions. We predict that two populations of B/PS bulges could exist and even coexist in the same galaxy. 



Key words. Galaxies: active 
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1. Introduction 
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How galactic bulges form is one of the leading questions for the 
gala xy formation theorie s. A consensus is begining to emerge 
(e.g. lAthanassoulal l2005h that two scenarios of bulge formation 
compete. The older scenario is the formation by an initial gravi- 
tational collapse or, in a more recent scenario by a series of mi- 
nor mergers in a similar way to massive elliptical galaxies. The 
second s cenario relies on the secular ev olution of stellar discs. 
This led iKormendv & Kennicutj (|2004|) to make a distinction 
between "pseudobulges" that are "bulges" formed through sec- 
ular evolution, and "classical" bulges, those with round smooth 
isophotes that show no discy structure in the central regions and 

thus are built up through mergers or a dissipative collapse. 

Th is dichotomy has been slightly refined by lAthanassoulal 
(1200 5*) who proposed to further split the pseudobulge category 
into two classes. The first one contains boxy-peanut shaped 
bulge (B/PS bulges hereafter) that are due to the vertical orbital 
structure of stel l ar bars seen edge -on ("Combes & Sande rslll98ll: 
|Pfennigei1ll984l:ICombes et al.|[l9 90; Pfenniger & Friedlilll99l|). 
The frequency of B/PS bulges is high: 45% of all bulges are 
B/PS, while amongst those the exact sh ape of the bulge depen ds 
mainly on the viewing angle to the bar dLiitticke et al Ll2000allbh . 
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The observed incidence of B/PS bulges is however consistent 
with that expected if they are as sociated with relativ ely strong 
bars. Near infrared observations (iBureau et al l l2006h have also 
revealed several B/PS features in 'classical' bulges. As shown 
by N-body simulations, true peanuts are bars seen side-on, i.e. 
with the major-axis of the bar roughly perpendicular to the line- 
of-sight. For less favourable viewing angles, the bulge/bar looks 
boxy, and if the bar is seen end-on it looks almost spherical. 
Stronger bars also lead to more prominent peanut shapes, as 
demonstrated ob servationally (e.g. [Lti tticke et al., 2000b) and 
theoretically (e.g. lBureau & Athanassoulall2005h . 

According to l Athanassoulal ( l2005h , "disk-like" bulges (DL 
bulges hereafter) belong to a second class. They are formed by 
star formation occuiTing in the gaseous inflow possibly driven by 
a stellar bar. Bulges formed according to this scenario can have 
observational properties attributed normally to stellar discs (ex- 
ponential photometric profiles, blue color, substructures like spi- 
ral arms, nuclear bars, circumnuclear rings, etc). In general, they 
can contain a measurable amount of gas, as well as a young stel- 
lar population sometimes distributed in bright spots. According 
to their mode of formation, DL bulges should have a much 
smaller scaleheight than B/PS or classical bulges since the gas 
distribution is rather flat even in the central galactic region. It 
is thus very q uestionable to name these structures "bulge" (see 
discussion b v lAthanassoulal r2005h . 
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^From the simulation point of view, it has been con- 
sidered that DL bulges deserved little attention. However, 
they did appear under several other names (exponential 
bulges, nuclear discs, double-bars, bar-within-bar, etc) in 
a number of numerical simulations or theoretical work s 
JShlosman & Begelman, 1989; Shlosman e t"an. ll99C 
Friedli & Martinet, 1993; Friedli et all 1 1996c iWozniak et al 



Table 1. List of runs. Masses are in 10' ^M0 units. 



2003; Wozniak & Champavert, 2006*). The main difficulty in 
self-consistently studying the formation mechanism of such 
bulges or other related central/nuclear morphological features 
comes from the fact that numerical codes have to include a 
dissipative component and recipes to mimic star formation and 
feedback processes, even using very simple rules and/or crude 
approximations. Apart from some peculiar cases, collisionless 
N-body codes are thus not able to reproduce DL bulges. 

DL bulges could be asso ciated with central velocity dis- 
persion drops (Wozniak et al., 2003) such as thos e observed by 
lEmsellem et al. (200D and Mai-quez et al.1 (|2003|) . The cr-drops 
result from the concentration of new stars toward the centre, and 
because this population of new stars is newly formed from the 
low-dispersion gas component, its velocity dispersion is much 
lower than for the old population. This effect is also amplified 
by the fact that the gas dispersion also drops toward the centre 
(and therefore the new stellar component too). This is due to the 
strong accumulation of gas toward the centre, in a nuclear disk, 
where dissipation is stronger than elsewhere, and where the gas 
therefore cools down efficiently. The stellar velocity dispersion 
could remain low even if the s tar formation rate is rather low 
jWozniak & Champaverlll2006l) . e.g. 1 Mq yr that is the same 
order of magnitude as the typical gaseous mass inflow rate into 
nuclear rings (Regan, Vogel & Teuben, 1997). This is additional 
evidence of a small DL bulge scaleheight. 

We decided to investigate the edge-on properties of B/PS and 
DL bulges. As stated above, this kind of study must perform 
N-body simulations including gas and star formation recipes. 
Typical simulations are presented in Sect. |2l In the rest of this 
paper we mainly discuss the formation mechanisms of both the 
B/PS and DL bulges and their dynamical properties. In partic- 
ular, we show that the young stellar population is able to form 
a B/PS bulge via a slightly different dynamical mechanism than 
classical B/PS bulge (Sect.O. The formation of this young B/PS 
bulge is accompanied by a flat nuclear disc (Sect.|4|. This could 
be interpreted through a linear and non-linear dynamical analy- 
sis of the families of resonant orbits (Sect.|5]). 

It should be stressed that young populations are much 
brighter for a few 10^ yr than old populations, so a DL bulge may 
also result from the light contrast, and not only from the mass 
distribution. In accompanying papers we will therefore address 
the issu e of detecting youn g B/PS bulges predicted by our simu- 
lations dMichel-DansSr&^Wozniakl |2008) and their kinematics 
(in preparation). For this purpose, we will use photometrically 
calibrated simulations including the absorption by the dust dis- 
tribution in the disc. 

2. Description of the numerical simulations 

For clarity, we will concentrate on a single case extracted from 
a dozen such simulations of various resolutions and initial se- 
tups. The generic simulation, named A^* hereafter, is thus rep- 
resentative of our database. Othe r examples can be found in 
iMichel-Dansac & WozniaS ( l2008h . 

An initial stellar population is set up to reproduce a typ- 
ical disc galaxy. Positions and velocities for 2.5 10^ parti- 
cles are drawn from a superposition of two axisymmetri- 
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cal 'Miva moto & Nagail (Il975h discs of mass Mi and M2 (cf. 
Table [T]), of scalelengths 1 and 3.5 kpc and a common scale- 
height of 0.5 kpc. Initial velocity dispersions are com puted solv- 
ing numerically the Jeans equations according to the Hernquis^ 
( 1993) method. The method was extended to take account of the 
presence of a dissipative component when solving for the stel- 
lar equation. The initial velocity dispersion was chosen to be 
anisotropic with (Tr = cr^ and cTg - cr^/c^/(40^), where cr,., erg 
and cr- are three components of the velocity dispersion along re- 
spectively the radial, azimutal and vertical directions and k and 
Q are respectively the radial and angular epicyclic frequencies. 
The resulting initial Q parameter in the central 500 pc radius 
increases quickly with radius from 1 to 1.5 and then slowly in- 
creases up to 2.3 in the external part of the disc. 

The initial disc radius is 30 kpc. The gaseous component of 
run A^^ is represented by 50000 particles for a total mass of 
1.1 10'*' Mq distributed in a 6 kpc scalelength Miyamoto-Nagai 
disc. 

A reference run of pure collisionless particles, named B"°^^, 
has been computed to carry out various dynamical comparisons. 
For this homologous run. Mi and M2 have been proportionally 
scaled so as to keep the same total mass and spatial distribution 
as A=^. 

The evolution is computed with a particle-mesh N-body 
code, derived from the original version of th e Gene va group 
(Pfenniger&FriedH, 1993; Friedli & Benj Il993h . which 
includes stars, gas and recipes to simulate star formation. The 
broad outline of the code is the following: the gravitational 
forces are computed with a particle-mesh method using a 3D 
log-polar grid with (Nr,N^,Nz) - (60,64,312) active cells. 
The smallest radial cell in the central region is 36 pc large 
and the vertical sampling is 50 pc. The extent of the mesh is 
100 kpc in radius and +7.8 kpc in height. The hydrodynamics 
equations are solved using the SPH technique. Since we used a 
polar grid and we need an accurate determination of the forces 
in the central region, we have improved the pre-computation of 
self-forces by subdividing each cell in {nr,n^,n^) = (32,6,6) 
subcells. Self-forces are then linearly interpolated before 
being subtracted from the gravitational forces. The spatial 
resolution and force accuracy are thus much high er than in any 



of ou r previous studies bas ed on the same co de dFriedli et al, 



119961: IWozniak et al.L l2C)03l: Michel-Dansac & Wozniak, 2004 
Hernandez et al.', '2005"; Michel-Dansac & Woznia^ i2006 
Emsellem et al., 2006; Wozniak & Champavert, 200^. 

The star formation process is based on Toomre's cri- 
terion for the radial instabili ty of gaseous discs (cf. 
iMichel-Dansac & Woznial |2004| for more details). When 
star formation is active, the radiative cooling of the gas is 
computed assuming a solar metallicity. In Figs. [T] and |2] we 
display the face-on and edge-on views of the bar region of A^^ 
for four times which will be used throughout this paper. They 
have been chosen as being illustrative of the bulge evolution. At 
the end of the simulation (f w 3000 Myr), the total number of 
particles is roughly 3.2 10^ for the stellar component and 30 000 
for the gaseous one. 45% of the gas has been transformed into 
stellar particles, mainly in the central 10 kpc. 
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Fig. 1. Face-on projected mass distribution in the central 10 kpc for t - 600, 1000, 1500 and 2000 Myr. The snapshots have been 
rotated to align the bar with the x-axis. 'young pop.' means only the population of particles created during the run. 



The main effect of a live dark halo (except to flatten the ro- 
tation curve of the disc at a large distance) is to permit the ex- 
change of angular momentum with the stellar disc. The rate and 
the amplitude of these exchanges depend on the velocity dis- 
persion of both the disc and the halo, and on the relative halo 
mass (e.g. Debattista & Sellwood, 2000; Athanassoula, 200^1 
FValenzuela & Klvpinll2003l) . The stellar disc could lose between 
a few % and 40% of its angular momentum mainly through reso- 
nances. Depending on the rate at which the stellar disc losses its 
angular momentum, the bar grows quite differently. Considering 
[Martinez- Valpuesta et al. (2006) simulations as representative, 
roughly 2/3 of the angular momentum loses by the bar-unstable 
part of the stellar disc is absorbed by the halo, the rest going to 
the outer disc. Most of these exchanges happen during the buck- 
ling of the bar. Afterward, the halo absorbs all the angular mo- 
mentum lost by the disc, l eading to a second phase of buckling 
(iMartinez- Valpuesta et ali l2006). The lack of a live dark halo in 
our simulations thus has the main consequence that we are not 
able to find any second buckling phase. 

However, our simulations are not completely devoid of ver- 
tical exchanges since the Miyamoto-Nagai density distribution 
allows us to build substantially inflated bulges, which is the case 
with our choice of parameters that leads to an SO-like initial stel- 
lar distribution. Thus the evolution of the particles confined near 
the z - plane is also driven by vertical exchange of angular 
momentum. Moreover, this paper is a first report of the eff'ect of 
a gaseous component with active star formation on the vertical 
structure of bars. It is important to be able to disentangle the ef- 
fect of the dissipative component from that of the halo on the 
disc evolution. Thus, simulations with a live dark halo, that are 
still in progress, will be reported elsewhere. 



3. Central region thickening 

In the case of pure N-body simulations, i.e. made of coUi- 
sionless particles only, whenever a disc galaxy forms a bar, 
a B/PS bulge develops in a few dynamical times. This pro- 



cess was firstly studied by Combes & Sanders! (11981 ). con 
firmed by ICombes et al] (I l990h and lRahaetal.1 (1199 lb . and 



recently reanalyzed by many authors (iMi hqs etali 1995 
Athanassoula & Misiriotis, 2002; Athanassoulal l2002l 12003 
i2005i ; iMartinez- Valpuesta et aU i2006 , etc). A stellar bar starts 
thin and then buckles out of the plane after a short phase of mir- 
ror symmetry breaking in the z direction (IPfenniger & Friedlil 
1991.) . Thus, initially asymmetrical with respect to the 
equatorial plane, the bar fina lly tends towards symmetry 
( IMartinez- Valpues ta et al.. 2006). Only a fraction of the whole 
bar buckles dAthan assoula, 2005), so that the outermost part of 
the bar stays thin and planar. 

This is once more what happens for the bar in B"°^^ (see 
Fig. [3]). Between t a 1400 and 1600 Myr the vertical thickening 
starts due to a break in the z mirror symmetry. It is only for t > 
1600 Myr that the z-distribution becomes more symmetrical. 

In the case of A^^ most of the young stellar population 
lies in a razor-like central disc during the first 450 Myr. This 
is due to the small vertical scaleheight of the gas distribution 
that remains thin because of its dissipative nature. It is well- 
known that such a stellar razor-thin disc is highly unstable (e.g. 
Merritt & Sellwood, 1994). Indeed, 450 Myr after the begin- 
ning of the young disc formation, the most central part of the 
disc starts to thicken out of the equatorial plane. In roughly a 
bar rotation period, the vertical distribution becomes symmetri- 
cally peanut shaped over the central 2 kpc (Fig.|2l right panel, at 
t - 600 Myr), while the young bar is approximately 8 kpc long 
(Fig.lTJ. However, at this time, the total mass of the central disc 
still being low, the thickening process has no detectable effect on 
the global mass distribution (Fig.|2] left panel). 

Afterward, the peanut-shape widens out as the young disc 
continuously evolves and increases in mass. The thick part of 
the disc doubles its radial size in less than 1 Gyr. The vertical 
scaleheight also increases with time leading to a well-developed 
peanut-shaped bulge for f > 1500 Myr. It is noteworthy that 
the total mass distribution, hence including both the initial and 
young population, plainly displays a B/PS bulge. This is mainly 
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Fig. 2. Left panel: edge-on projected mass distribution in the central 20x6 kpc for t - 600, 1000, 1500 and 2000 Myr The snapshots 
have been rotated to align the bar with the x-axis. Right panel: only the population of particles created during the run. 



due to the mass of the young population that amounts to a sig- 
nificant fraction of the whole mass at that time. Indeed, the B/PS 
bulge is clearly much less marked if we only display the mass 
distribution of the initial population (Fig.|4]i. 

Disentangling the two causes that could be responsible for 
the increase in boxiness of the initial population is not so ob- 
vious. Indeed, as the mass of the young population trapped in 
the B/PS increases, the gravitational potential deviates more and 
more from its initial quasi-spherical shape. This leads to a ver- 
tical mass redistribution of the initial population since particles 
can be trapped by orbits associated with vertical resonant fam- 
ilies. But the initial disc p opulation could b e also unstable to- 
wards vertical instabilities (Raha et al., 1991). This is supported 
by looking at the homologous simulations B"°^* that also develop 
a B/PS bulge (Fig. [3]). For B"°^^ we said above that the growth 
mechanism is the same as for all other collisionless simulations 
of cold stellar discs performed so far 

However, a few differences between B"°^^ and A^* should be 
noticed: 

1 . At the end of the buckling phase, the extent of the boxy re- 
gion is slighdy wider for B"°^* than A^^. Moreover, B"°^^ 
is much more peanut-shaped than A^^, especially at f = 
2000 Myr and beyond. In Sect. |5] we attempt to explain this 
in terms of resonances. 

2. The buckling is clearly asymmetrical for B"°^^ between f - 
1400 and 1600 Myr, whereas it remains symmetrical for A^^. 
If any asymmetry in the mass distribution appears during 
the evolution, its scaleheight should be less than the ver- 



tical resolution of the code (i.e. 50 pc). Indeed, since the 
z-distribution of the A^^ young population remains always 
symmetrical with respect to the equatorial plane, the poten- 
tial well created by the young population is also permanently 
symmetrical. It thus could be more difficult to break the z- 
symmetry for A^^ than for B"°^^. This difference supports the 
interpretation that vertical resonant orbits are populated. 



In principle, we cannot exclude the possibility that the thick- 
ening of the young population might be asymmetrical because 
low-order bending modes, those that can give an asymmetrical- 
shaped mass distribution, can be eras ed by our 50 pc v ertical res- 
olution. Indeed, it has been shown bv lMerritt & Sellw ood ( 199^) 
that finite grid effects, especially the vertical mesh resolution, 
could stabilize the low-order bending modes in the case of a per- 
fectly thin disc. However, the young disc created in all our simu- 
lations is always embedded in the initial stellar population which 
has a much higher scaleheight (initially 500 pc, increasing with 
time). The vertical distribution of the initial stellar population 
also remains symmetrical with respect to the equatorial plane 
(cf. Figs. [3] and |4]i. Moreover our vertical resolution is much 
higher than for the Merritt & Sellwood simulations. The verti- 
cal oscillating frequencies of particl es thus are not domin ated by 
the mesh spacing as is the case for lMerritt & SellwoodI (119941) 
simulations. 
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Fig. 3. Edge-on projected mass distribution of run B"°^^ in the 
central 20x6 kpc for t = 600, 1000, 1500 and 2000 Myr. 
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Fig. 4. Edge-on projected mass distribution of the initial pop- 
ulation in the central 20x6 kpc for t =600, 1000, 1500 and 
2000 Myr. 
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Fig. 5. Face-on gas mass distribution in the central region where 
the flat nuclear stellar disc made of the young population devel- 
ops for A^^. The gas particle distribution has not been convolved 
by the SPH kernel to emphasize the ring structure. 

4. Flat nuclear disc 

As discussed above, a thin stellar disc cannot remain flat for a 
long time, such as the A^^ initial young central disc. However, for 
A^^ even if a significant part of the young stellar disc thickens 
rapidly, there is a small nuclear component which remains flat 
for a very long time (Fig.lSj, i.e. until the end of the simulation 
3.5 Gyr). This nuclear component is made up of a young stellar 
population formed in the gaseous disc accumulated in that region 
because of the torques exerted by the large-scale stellar bar. The 
formation of this flat nuclear structure is thus concomitant with 
the thickening of the central part discussed in Sect. |3] After 2 
Gyr, the stellar mass amounts to roughly 7.25 x 10^ Mq in a 
cylinder of R - 0.5 kpc and z + 0.5 kpc. The young population 
accounts for roughly 34% of the total mass. 

The nuclear disc quickly develops a small bar with its own 
pattern speed, almost 10 times higher than the large-scale bar 
The nuclear bar is encircled by a circumnuclear ring (Fig. |5]l. 
The nuclear bar appears at f « 575 Myr and is long-lived, al- 
though at the resolution of our simulations (between 36 pc at 
the centre and 100 pc at 1 kpc), it periodically dissolves in a 
spiral-like structure. The detailed study of the morphological, 
dynamical and kinematical properties of the nuclear bar and the 
circumnuclear ring deserves a ded icated paper and thus will be 
publish elsewhere dWozniald l2008i) . Hereafter we ignore the in- 
ternal structure of the nuclear disc. 

The existence of flattened and rapidly rotating nuclear stel- 
lar discs has been pr edicted by [Sh losman &Begelman (1981) 
and further studied bv lShlosman et al.l (11990 ). They showed that 
such stellar discs could remain flattene d for a long time since the 
two-body relaxation is a slow process. [Chung & Bureaul (|2004|) 
detected such a nuclear disc in roughly 1/3 of their sample of 
24 edge-on B/PS galaxies. They have given clear evidence that 
these discs exhibit a cr-drop in their centre. They also speculated 
on the possibility that such discs are formed through gas inflow 
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Table 2. Rough sizes of the box (from rotation axis to corner) 
and the central stellar disc (radius in kpc) determined on the 
mass distribution of the young population alone. The radius of 
main linear resonances for A^^ and B"°^^. hILR is the horizontal 
inner Lindblad resonance, vILR the vertical one and hUHR is 
the horizontal ultra harmonic resonance. 
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Myr 


Box and disc sizes 










disc 


0.27 


0.43 


0.50 


0.50 


boxy 


0.67 


1.17 


1.40 


1.59 


boxy/disc 


2.48 


2.69 


2.80 


3.18 


Resonance radii of A^* 










hILR 


1.27 


1.61 


1.89 


1.80 


vILR 


N/A 


1.34 


1.68 


1.81 


liUHR 


3.57 


3.82 


4.21 


4.00 


Resonance radii of B"°^* 










liILR 


1.61 


2.24 


2.45 


2.84 


vILR 


N/A 


N/A 


N/A 


2.95 


liUHR 


4.69 


5.12 


5.25 


5.26 



and subsequent star formation. The whole initial young central 
disc discussed in Sect. [3] cannot be analogous to the flat disc of 
[Shl osman & Begelman ( 1989), but its nuclear component obvi- 
ously has similar properties. 

^From the kinematical point of view, the nuclear compo- 
nent is undou btedly associated with the cr-drop phenomenon. 
IWozniak et al.l (|2003) showed that such a young nuclear stel- 
lar disc is responsible for the cr-drop since it is the kinemati- 
cal signature of stars that have been born from a dynamically 
cold gaseous component. The presence of the nuclear bar and 
transient spiral structure in A^^ marginally increases the radial 
velocity dispersion, but the effect on the line-of-sight velocity 
dispersion remains weak so that a cr-drop should be visible in 
double-barred galaxies. A^^ thus confirms the potential relation- 
ship bet ween nuclear disc and cr -drop in B/PS bulges as sug- 
gested bv lChung & Bureaul (|2004 . 

The respective sizes of the nuclear disc and peanut-shaped 
bulge are given in Table |2] The two components enlarge as the 
galaxy evolves but the ratio of the box over the disc sizes also 
slightly increases with time. The two phenomena, whose the 
common initial cause is the presence of a young stellar popula- 
tion formed in the inflowing gaseous material driven by a large- 
scale bar, seem dynamically distinct. Both morphological struc- 
tures seem to evolve independently from each other However, 
the main driver of the internal dynamics is well-known to be the 
large-scale bar and the set of resonances associated with its ro- 
tation pattern. 

5. Horizontal and vertical resonances 

5.1. Linear analysis 

Numerous authors have tried to correlate the size of mor- 
phological structures to the dynamical resonance locations. 
Circumnuclear and outer rings seem to be correlated with the 
location of, respectively, the inner Lindbl ad reso nance (ILR) 
and outer Lindblad resonance (OLR) (cf. lButa & Combesll996h . 
The ratio of nuclear bar length to that of the large-scale bar could 
be similar to the ILR to corotation (CR) ratio, the nuclear bar 
corotation being dynamically coupled to the large-scale bar ILR 
(e.g. Rautiainen & Salo 1999). H owever, some other simulations 
did not show such coupling (e.g. lHeller et al.ll200lb : this matter 
is still under debate. 




Fig. 6. Resonance diagram (Q - Q.p)/K (horizontal resonances, 
full line) and (Q. - Op)/v (vertical resonances, dashed line) at 
f = 1500 Myr. Location of ILR (dot-dash), UHR (triple dot- 
dash) and CR (dot) are plotted. The shaded box depicts the re- 
gion occupied by the nuclear disc. 



To determined the location of the linear CR and ILR dy- 
namic resonances, we computed the circular orbit frequency 
and the radial epicyclic frequency k as dPfennigeit 1 1 990l) : 



^(^) - + — +2 {-—\ 

dx^ dy^ \r dr I 



(1) 



(2) 



where O is the gravitational potential and (• ■ ■) stands for an az- 
imuthal average. 

Strictly speaking, these frequencies predict the oscillation 
frequencies of the orbits in the axisymmetrical case only. They 
do not provide any indication of whether families of periodic or- 
bits do follow such oscillations when the bar growth breaks the 
axi symmetry. However, a number of previous orbital studies (cf . 
Mic hel-Dansac & Wozniakl l2006t IWozniak & Michel-Dansad, 
2007, and discussion therein) suggest that the epicyclic approx- 
imation could lead to an acceptable estimation of the resonance 
locations displayed in Table |2l in particular if we are mainly 
interested in their evolution rather than their accurate absolute 
position. For instance, usin g a careful integration of orbits to 
accurately compute Q. and k, Michel-Dansac & Wozniakl (l2006h 
found that the error on Rqr remains within 10%. This technique 
has been used to study a few snaphots when looking for higher 
order resonances (cf. Sect. 15. 2t . 

The flat nuclear disc is entirely inside the linear ILR of the 
large-scale bar, its radius being roughly 0.5 kpc at f = 1500 Myr 
Indeed, 7?hiLR the horizontal ILR (hILR) radius computed using 
the linear approximation is ^ 1 .9 kpc from the centre (cf. Fig.|6]l. 
The gaseous component which also forms a nuclear disc occu- 
pies the same region as the young stellar nuclear disc. During 
the evolution of the disc, it is well known that the large-scale 
bar slows down, leading 7?hiLR to increase with time. Hence one 
could imagine that the size of the nuclear disc increases propor- 
tionally to /?hiLR but in fact its growth saturates after roughly 
1 .5 Gyr. The nuclear disc does not entirely fill the region encir- 
cled by the hILR. Indeed, the radius where Q.- Kjlis maximum 
is expected to be the limiting dynamical radius. In Fig. |6] the 
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Fig. 7. Evolution of horizontal (lines) and vertical (dots) linear 
ILR radius for run A^^ (black thick line and full symbols) and 
gnosf ((Jotted lines and opened symbols). 



limiting radius is roughly where (Q. - Q.p)/Kis maximum, that is 
quite close to an Q - /(■/2 local maximum. 

The combined effects of the horizontal and vertical ILRs 
cause the planar orbits lyin g in the equatorial plane to be destabi- 
lized dCombes et al.L[T990b . but the B/PS re gion extends outside 
the IL Rs. However, as already outlined bv IPfenni ger & Friedlil 
(1199 Ih . in order to observe a well defined boundary in the 
B/PS it is necessary that families of orbits supporting the shape 
cease to exist or, at least, become unstab le beyond some well- 
defined height. iPfenniger & Friedlil ( Il99lh suggested by default 
that corotation could be the limiting resonance. However, at 
t = 1500 Myr, the vertical ILR (vILR) is located at « 1.7 kpc, 
which is very similar to the distance of the corner of the B/PS 
from the axis of rotation (cf. Table |2]i that could be the ex- 
tent of orbit families associated with the vILR. This obviously 
suggests that the ILRs, and their associated families of reso- 
nant orbits, are responsible for limiting the extent of orbits in 
the vertical direction. Since it has been shown that the horizon- 
tal UHR delineates the bulk of the bar in the equatorial plane 
(Michel- Dansac & Wozniakl 120061: IWozniak & Michel-Dansad 
12007 ). the combined effects of UHR and ILRs play a very im- 
portant role in shaping the central region of disc galaxies. 

At t = 1500 Myr (Fig. |6ll the location of the hILR is 
very close to the location of the vILR (Rwilr~ 1-7 kpc at f = 
1500 My r). This coincidence has been observed by pr evious au- 
thors (Combe s et an.ll990l:|Pfenn iger & Friedlil ll99lh and it has 
been suggested that this situation could have some consequences 
on the local dynamics. As shown in Fig. |7] for A^^ the coin- 
cidence occurs during all the run, /?v1lr being always slightly 
smaller than /?hiLR, although the difference decreases with time. 
The rapid changes in the nuclear mass distribution due to the gas 
inflow imply rapid fluctuations of the ILR positions as well, the 
fluctuation timescale typically being the local dynamical time 
(less than 50 Myr). 

For B"°^*, the coincidence arises progressively (cf. Fig. |7]). 
Indeed, during the first 1.6 Gyr there is no vILR. The vILR 
appears at a radius ,RvIlr~ 1-4 kpc and then increases to coin- 
cide with the hILR at f x 2.1 Gyr. We have said above that the 
B/PS bulge formation for B"°^^ follows a different path than for 
A^^. Indeed, before the vILR appears, the stellar disc buckles out 
asymmetrically. It is likely only when families of orbits associ- 
ated with the vILR appear that the mass distribution could be 



symmetrized with respect to the equatorial plane. This is also 
the case for A^^ although both ILRs appear simultaneously. The 
thickening of A^* is thus symmetric because ILR resonant fami- 
lies exist from the beginning of the bulge formation. 

5.2. Non-linear analysis 

To be able to discuss the potential effect of higher order reso- 
nances, we have computed the orbital frequencies Q, k and v of 
a representative sample of particles for a few selec ted snapshots. 
We ha ve applied a variant of the technique of lAthanassoulal 
(|2003|) . We have frozen the potential at a given time in the 
simulation and then computed orbits in the inertial potential 
to determine the principal f requencies using the technique of 
ICarpintero & Aguilail ( Il998l) . We used almost 50000 particles 
as initial conditions chosen at random in a limited domain of the 
phase-space. Indeed, being mainly focused on the resonances 
inside the bar, we restrict our computation to the particles that 
a priori reside in the region encircled by the corotation. We a 
posteriori checked that the (Ej, r) space (Ej being the Jacobi 
constant) and time-averaged Lindblad diagram (E j, L,) are well 
sampled by t his choice. To discuss our re sults we will used the 
notation of iSellwood & WiUdnsonl d 19931) for closed (periodic) 
orbits. In this notation, m:n:l implies m radial oscillations in the 
(x, y) plane and n vertical oscillations in z as the orbit achieves I 
rotations about the center. 

We have selected the snapshot t - 2000 Myr for B"°^* which 
is the moment when the B/PS shape is well established and the 
snapshot t - 1500 for A^*. To easily identify periodic orbits that 
shape the morphology of barred galaxies, we display the distri- 
bution of -TC = (Q - Q.i,)/k and // = (□- Qp)/v in Figs. H to 
|9] Families of periodic orbits thus have commensurable ratios 
7C = l/m and N = l/n. 

In these Figures, we have discarded orbits for which the max- 
imal absolute radial extent is less than 1 kpc to avoid contamina- 
tion by quasi-circular orbits that remain confined in the central 
nuclear disc. These flat orbits do not participate in shaping the 
bulge and mostly populate the regions 'K = I and N - 0. On the 
contrary we did not try to filter out multi-periodic orbits (/ > 1) 
so that they are mixed with I - 1 orbits. 

For the two snapshots, the dominant family of orbits is re- 
lated to the equatorial 2:n:l resonance, that is the external hILR 
(TC = 1 /2) since we have discarded as much as possible resonant 
families associated with the innermost hILR. However, both the 
'K and the TV diagrams are quite different for A^* and B"°^^ and 
deserve to be separately discussed, then compared. 

For B"°^^ (Fig.ID, the dominant family is made of 2:2:1 or- 
bits (TC = 1/2 and N - 1/2). These 2:2:1 or bits are known 
to be responsible for the classic al B/PS bulges dCombes et al.L 
1 1 9901: IPfenniger & Friedlil 1 1 99 Ih found in coUisionless N-body 
simulations and ha ve thus focused m ost att ention under variou s 
names (e. g. B^, for 'Pfenniger'(1984'), zl for lHasan etan ( ll993h . 
xlv3 for 'Skokos et al. (2002)). These orbits dominate the 7C 
and N distributions. However, since the mass fraction trapped 
around the horizontal and the vertical frequencies are different, 
this suggests that a fraction of the TC - 0.5 peak is populated by 
2:«:1 orbits with n + 2. 

Several other non-negligible contributions are concentrated 
around -TC = 0.25 (4:«:1 orbits related to the UHR), 0.4, * 0.6 
and 0.75. It is much more difficult to uniquely identify the 
orbit families responsible for these peaks but they mostly contain 
multi-periodic trapped orbits. For instance we found 5:n:2 orbits 
that contribute to the = 0.4 peak. 
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Fig. 8. Distribution of mass as a function of '7C = (Q - Q.p)/K 
(upper panel) and N - (Q. - ilp)/v (lower panel) for B"°^^ at 
t - 2000 Myr The total mass is normalised to unity. The 2:1:1 
CTC = 0.5), 3:1:1 (-TC = 0.33) and 4: 1:1 (TC = 0.25) resonances are 
plotted as dotted lines. Multi-periodic orbits (Z > 1) are super- 
posed on / = 1 orbits. 
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Fig. 9. As Fig. IDfor at t = 1500. The peak at TC = 0.5 has 
been truncated to 0.15 but extends to a value of 0.20. 



For A^^ (Fig. |9l), the distributions are very different even if 
the main contribution is still concentrated around "K - 0.5 as for 
gnosf jjjg fraction of mass (roughly 0.2) trapped for TC = 0.5 
is however much higher than for B"°^^ (0.13). The secondary 
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Fig. 10. Examples of dominant families of orbits for at 
t=1500 Myr. From top to bottom, a quasi 2:2:1 orbit and three 
different instances of 2:3:1 orbits. 



peaks are less prominent than for B"°^^ and they form a forest of 
frequencies betwen 'K = 0.1 and 0.8. 

The vertical structure is also very different from the colli- 
sionless run. The mass trapped around N = 1/3 (the vertical 3:1 
resonance) is slightly greater than for N = 1/2 resonant orbits. 
TC = 1/2 for these orbits, as for 2:2:1 orbits, so that they must to 
be classified as 2:3:1 resonant orbits. A few examples are shown 
inFig.[IO] 

The role of this family of orbits has been emphasized by 
iHeller & Shlosmanl (Il996h who demonstrated that it is con- 
nected to and amplified by the p resence of a massive nuclear 
ring. Indeed, Hel ler & Shlosmanl ([1996) showed that the ma- 
jor periodic orbits within the corotation radius are affected by 
the perturbation of a massive ring. The main family of or- 
bits that sustain a stellar bar (the so-c alled x\ family in the 
IContopoulos & Papavannopoulosl ( ll980l) notation) becomes ver- 
tically unstable for several ranges of Jacobi constant (E j) values. 
In the E j domain where xi coexists with X2 and X3 families (that 
are respectively stable and unstable elliptical-like orbits perpen- 
dicular to the bar and that occupy the region between the two 
ILRs), two xi instability strips provoke z and z bifurcations to 
2:2:1 families (the so-called symmetrical BAN - banana - and 
asymmetrical ABAN - anti-banana - families) and 2:3:1 fami- 
lies^ 

iHeller & Shlosmanl (119961) showed that the width of these in- 
stability gaps depends on the mass of the ring. As the ring mass is 
increased, the instability strip responsible for 2:3:1 orbits grows 
in size and moves toward lower Ej. In some extreme situations 
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Fig. 11. Examples of other boxy-shaped famiUes of orbits for 
at t=1500 Myr. 



where the mass of the ring is very high (e.g. 1O^M0) some in- 
stabihty strips appear for the X2 family. 

2:3:1 orbits are symmetrical with respect to the (x,z) 
plane and antisymmetrical with respect to the (y,z) one (or 
the converse) unlike 2:2:1 orbits that are either symmetri- 
cal or antisymmetr i cal ab out both the ix,y) and (y,z) planes. 
[Heller & Shlosman ( 1996) (their Fig. 7 and 8) and Skokos et alj 
(2002) (their Fig. 9) display typ i cal m embe rs of this family. I t 
has been noted bv lSkokos etalj (|2002|) and lPatsis et"!!] (l2002h 
that the 2:3:1 family (called Ji:lv3 in these papers) as well as other 
families associated with higher order vertical resonances (i.e. 
2:4:1, 2:5:1 etc.), could play an important role in shaping small 
sized B/PS bulges because the z extent remains small whatever 
the values of Ej are. For A^^ they appear to be a major contribu- 
tion to the shape of the young B/PS bulge whereas they are fully 
absent for B"°=^. 

For A^* the main family after 2:2:1 and 2:3:1 i s the 2:5:1, 
visible in Fig. |9] around N — 0.2. As discussed bv iPatsis et alj 
(I2OO2I) this family is also a strong contributor to the B/PS bulge. 
But, as for 2:3:1 orbits, they are absent in B"°^^. Many other 
kinds of resonant orbits contribute to the B/PS. In Fig.[TT]we dis- 
play a few examples. Their contribution is mainly concentrated 
around TC ~ 0.2 and 0.66 whereas N - 0.2. 



Even if the famihes of closed periodic orbits are the back- 
bone of any stellar bar, strictly speaking they occupy a null vol- 
ume of the phase-space. Only orbits trapped around the stable 
families are responsible for the shape of the bar. These trapped 
orbits have fundamental frequencies slightly shifted from those 
of their parent families so that they are responsible for the broad- 
ening of the spectral lines around the commensurable value of "TC 
and N. 

6. Conclusions 

We have numerically investigated some dynamical properties 
of edge-on boxy or peanut-shaped (B/PS) and disc-like (DL) 
(pseudo-)bulges. Our results are summarised as follows: 

1. We are able to confirm that the B/PS in N-body collision- 
less simulations is due to the classical break in the z mir- 
ror symmetry. However, in our numerical simulation that in- 
cludes a gaseous component and star formation recipes, the 
bulge-growing mechanism is quite different from the pure 
N-body case. The young stellar population that is born in a 
thin gaseous disc rapidly populates vertical resonant orbits 
triggered by the combined effects of the horizontal and ver- 
tical ILRs. This leads to a B/PS bulge mainly made of stellar 
material younger than the surrounding population. The mor- 
phology and extent of young B/PS bulges are significantly 
different from the classical B/PS bulge. We thus predict that 
two populations of B/PS bulges could exist and even coexist. 
They might be distinguished by deep photometric observa- 
tions or careful stellar population analyses. 

2. In N-body collisionless simulations the main orbit family re- 
sponsible for the B/PS is the 2:2: 1 . On the contrary, if a dis- 
sipative component is present and can form new stars, addi- 
tional asymmetrical families contribute to the B/PS. In the 
case of our simulation, 2:3:1 and 2:5:1 orbits trap a signifi- 
cant fraction of the mass. Their appearance could be linked 
to the massive circumnuclear ring. 

3. Aflat discy stellar component appears simultaneously with 
the thickening of the young population. It is due to star for- 
mation in the nuclear gaseous disc. Remarkably, it remains 
flat throughout the sim ulation (3.5 Gyr) although it de velops 
a bar, as predicted bv IShlosman & BegelmanI (Il989h . This 
suggests that nuclear bars, embedded in large-scale coun- 
terparts, are flat structures. Both the stellar and gas discs 
are located well inside the ILR, limited by the radius where 
(Q - Q.p)/K is maximum. 
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